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A REDUCED ORDER MODEL OF THE LINEARIZED INCOMPRESSIBLE NAVIER-STOKES 
EQUATIONS FOR THE SENSOR/ACTUATOR PLACEMENT PROBLEM 

BRIAN G. ALLAN * 

Abstract. A reduced order modeling approach of the Navier-Stokes equations is presented for the design of a 
distributed optimal feedback kernel. This approach is based on a Krylov subspace method where significant modes 
of the flow are captured in the model. This model is then used in an optimal feedback control design where sensing 
and actuation is performed on the entire flow field. This control design approach yields an optimal feedback kernel 
which provides insight into the placement of sensors and actuators in the flow field. As an evaluation of this approach, 
a two-dimensional shear layer and driven cavity flow are investigated. 

Key words. Navier-Stokes equations, feedback control, sensor/actuator placement, active flow control 
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1. Introduction. Experimental results using micro actuators and sensors have demonstrated that active flow 
control has the potential to increase the performance of high-lift wings, cavity noise, and other flow systems [13, 1 8, 5). 
Placement of these micro actuators and sensors can have a dramatic effect on the performance of the flow control 
system. The location of the sensors and actuators is complicated by the distributed nature of the nonlinear flow 
equations and the dynamics of the closed-loop system. The goal of this investigation is to identify the spatial regions 
of the flow where sensing and actuation are favorable under feedback control. 

By assuming control and sensing everywhere in the flow field, a distributed optimal feedback kernel can be 
computed. Evaluation of this feedback kernel shows spatial regions of the flow field which are more significant, 
in terms of actuation and sensing, than other regions of the flow field 1 1 1|. By identifying these regions of the 
flow, the search space for the placement of sensors and actuators can be reduced. This methodology does not give 
exact locations for point actuators/sensors and should be considered as a prefilter to a point actuator/sensor location 
placement problem. 

Calculation of the optimal feedback gains requires a finite dimensional approximation to the infinite dimensional 
Riccati equations. The cost of computing this finite dimensional Riccati solution is order N* where N is the number of 
states in the finite dimensional approximation to the dynamical system. Since the number of states needed to resolve 
a flow field can become very large, the cost of computing a solution to the Riccati equations can grow prohibitively 
expensive. However, this cost can be reduced by developing a reduced order model which contains significant dynam- 
ics of the flow system. To generate this reduced order model, the linearized incompressible Navier-Stokes equations 
are projected onto a Krylov subspace. This reduced order model is then incorporated into an optimal feedback control 
design. To evaluate this methodology a two-dimensional shear flow and a driven cavity flow problem are investigated. 

2. Governing Equations. This section describes the equations governing the dynamics of a viscous flow in 
two-dimensions. In this investigation the incompressible Navier-Stokes equation are represented in vorticity stream 
function form. This form of the Navier-Stokes equation was chosen over the primitive variable form, i.e. velocity and 
pressure, because it has a reduced number of unknowns with no incompressibility condition. A disadvantage to the 
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vorticity stream function form is that the no-slip boundary conditions produce two boundary conditions on the stream 
function and none for vorticity [16]. 

2.1. Two-Dimensional Navier-Stokes Equations. This investigation considers a two-dimensional square do- 
main ‘D with a boundary d‘D. The governing equations are the incompressible Navier-Stokes equations. These 
equations can be expressed in a vorticity stream function form. This form is achieved by taking the curl of the 
nondimensional momentum equation and substituting the definition of vorticity and stream function. The vorticity 
stream function form of the incompressible Navier-Stokes equations is then expressed as 

dco . 1 ^ 

— = — (u-V)oo-i--— V^to in ‘D 
dr Re 

V 2 \|/ = -0) in (D 


u = Ub on d*D 

where 0), \| /, and u are the dimensionless scalar vorticity, stream function, and velocity vector respectively. Note that 
the velocity at the boundary must satisfy a compatibility condition which follows from the integration of the continuity 
equation over the domain ‘D. Integrating the continuity equation over (D , applying the divergence theorem, and using 
the velocity boundary condition, results in 

(2.2) [ V u= <f Ub-n = Q 

J'l) Jd'D 

where n is the unit vector normal to d'D [16). The scalar vorticity field 0) is the z-component of the vorticity vector 
to— V x u which is given by 

(2.3) 0) = e- • (V x u) 


where the velocity vector u(x,f) — (w,v) and e- is the unit vector normal to the xy plane. The stream function \|/ is 
defined such that 


(2.4) 



v = 


3 \|/ 

a7 


This definition of the stream function produces a velocity field which exactly satisfies the incompressibility condition 
V • u = 0 for two-dimensional flows. The relation between the stream function and velocity vector can also be written 
in the compact form u = Vig x e z . 


2.2. Linearized Navier-Stokes Equations. The optimal feedback control design in this investigation is base on 
linear quadratic regulator (LQR) theory. Since the LQR control design is based on classical linear dynamic system 
theory, the nonlinear governing equations in Eq. (2. 1 ) can not be used for the feedback control design. Therefore the 
equations in Eq. (2. 1 ) are linearized about a desired base flow state where it is assumed that the flow is stabilized by the 
controller resulting in small perturbations about the base flow. These linearized Navier-Stokes equations describe the 
linear evolution of small perturbations about a given base flow field. It is desired that this base flow field be a steady 
state solution to the nonlinear governing equations. By making the base flow a steady state solution to the nonlinear 
equations, the time derivative of the base flow will drop out when the equations are linearized. 

This perturbation of the flow variables about some base flow state can be expressed as 

o)(x,/) = Q(x) +ea)(x,r) 


(2.5) 


<j/(x,0 = 4'(x) + e\|/(x,r) 


u(x,r) = U(x) + eu(x,/) 



where the flow states <5, vj>, and u satisfy equations in Eq. (2. 1 ) and e is some 'small' parameter. The perturbed states 
are represented by the variables to. \ I. and u and the base flow by Q. M*. and U = ( U,V ). Here the base flow states are 
assumed to he steady state solutions and are not function of time. 

Substituting Eq. (2.5) into Eq. (2. 1 ) and considering the terms which are of 0(e°) results in the equations, 

0= -(U-V)£2 + -|-V 2 Q in D 
Re 

(2.6) V 2 V = -Q in 2) 


U = V h on 3D 

These equations are the nonlinear steady state incompressible viscous flow equations. The desired flow state which is 
to be stabilized by the optimal feedback controller must satisfy the equations in Eq. (2.6). Considering now only the 
terms which are of 0(e) gives the equations, 

^ = — (U-V)co— (u V)£2 + -J-V 2 o) in© 
of Re 

(2.7) V 2 \|/=— to in (D 

u = uy, on d(D 

These equations describe the linear evolution of the perturbed flow stale about the base flow state. The nonlinear term 
(u ■ V)oo is of 0(e 2 ) and is dropped from above equations. 

The equations in Eq. (2.7) can be rewritten in the following conservative form, 

^ = -V-(coU)- V-(£2u) + -J-V 2 o> in© 
at Re 

(2-8) V 2 i|/ = -co in D 


u = u/, on d(D 

where the fact that the base flow velocity and the perturbed flow velocity satisfy the continuity equations V • U = 0 and 
V ■ u = 0 has been used. 

2.3. Boundary Conditions. The physical boundary conditions on the velocity u at a wall result in two boundary 
conditions on the stream function V|/. These two boundary conditions are derived by separating the normal and tangen- 
tial components of the velocity u at the boundary. Quartapelle (1993) shows that the boundary condition u| f >/> = u /, 
results in the two boundary conditions 

(2.9a) = a 


(2.9b) 



= b 


where a(s,f) = n(^) • u and b = -x • u/,. The variable .v is a coordinate along the boundary and si is 

any fixed point along the boundary. The vector x is a unit vector tangential to the boundary. For the flows considered 
in this investigation, the variable a(s,t) = 0. 

The two boundary conditions used in this investigation are periodic and no-slip nonporous walls. On the no-slip 
walls the Dirichlet boundary condition in Eq. (2.9a) was used to solve the Poisson equation for the stream function. 
The Neuman boundary condition Eq. (2.9b) for the stream function was used to solve the vortieity transport equation 
by deriving a Dirichlet boundary condition for vortieity. 
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3. Control Problem. The classical theory of control systems was developed for systems governed by finite 
dimensional ordinary differential equations (ODE), also known as lumped parameter systems. In this investigation the 
governing equations for fluid dynamics are partial differential equations (PDE) where the state of the system lies in 
some infinite dimensional function space. A system of this type is also known as a distributed parameter system where 
the states and control inputs are distributed spatially [10]. This section describes the LQR feedback control problem 
for the infinite dimensional system and its finite dimensional approximation. 


3.1. Infinite Dimensional Problem. The distributed control applied to the flow takes the form of a spatially 
distributed body force f(x,r). The torque generated by the body force, normal to the xy plane, is given by g — 
(V x f) *e-. This distributed control torque g(xj) appears on the right-hand side of the linearized vorticiiy transport 
equation. 


5(0 


-Nu to - Nqu + 


Re 


L(o + g 


in £> 


(3.1) 


via® = 


L\\f = —to in iD 



on d*D 


These PDEs are now written in an abstract form which is conducive to the optimal feedback control design. Treating 
this distributed parameter system as an evolution ODE, Eq. (3. 1 ) can be written as 


(3.2) 


ti = A($+Bg 


Cl)(x,0) = (Do 


where (i) = and A and B are infinite dimensional operators. The operator# for the system given in Eq. (3. 1 ) 

is just the identity but is included for completeness. Given a flow field (0, the action of the operator A on co is 

(3.3) = -V • (coU) - V • (ilu) + — V 2 w 

Re 

where the velocity u is an explicit function of the stream function vj/, which is an implicit function of vorticity co. 

The control input g is computed by using the state to(x,/) in the following way 

(3.4) g(xj)= [ K(x£)w&t)d$ 

J‘D 

where is the distributed feedback kernel. Therefore the control input g, at a position x, is given by an integral 

over the entire domain of the distributed gain K multiplied by the current vorticity field. 

The goal of this control design is to find a feedback kernel K which produced an optimal control input g opt . In 
order to define what optimal means, consider the following quadratic equation 

(3.5) j{a a ,g)= r[(<o(^0,e^)ft)(^0) + {^0,^)g(^0)]^ 

Jo 

where (7(x) and R(x) are weighting functions and the notation < *, ■ > represents an inner product over the domain T>. 
The term optimal is now defined as the control input which minimizes the quadratic cost function J . By minimizing 
7, the perturbations in the vorticity from the desired state £1 are minimize and the needed control input g is also 
minimized. 

The LQR problem for the infinite dimensional system is stated as 


(3.6) 


min7(o}(x,0),g(x,Q) 

s (*,0 
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subject to the system in Eq. (3.2). The optimal control input g opt (\j) which minimizes Eq. (3.6) is given by 

(3.7) gopi(x,t) = [ (0(^,0^ 

J'D 

where B* is the adjoint of B and FI is the nonnegative self-adjoint solution to the steady state, infinite dimensional, 
algebraic Riccati equation (ARE) 

(3.8) A'n + nA-nfl/r'/rn + c^o 

A solution n exists for Eq. (3.8) if the pair (A,B) is stabiliz.able and the pair (A,C) is detectable [4]. These two 
conditions are satisfied since actuation and sensing are performed everywhere in the flow field. From Eq. (3.7) it can 
be seen that the optimal feedback kernel K opt is given by 

(3.9) K„ rl =-R-'B* n 

3.2. Finite Dimensional Approximation. A finite dimensional approximation to the infinite dimensional feed- 
back kernel K is made using a finite difference method. In this approximation it is desired that the feedback kernel 
— y K, in an appropriate sense, as N — > The computational issues associated with this problem have been 

addressed by Banks et al. [2] and Gibson [8|. 

The finite dimensional approximation to Eq. (3.5) is expressed as 

c o N = A n to ;V + Bg N 

(3.10) 

0 /( 0 ) = < 

where o/ is a vector containing the spatially discrete values of the vorticity field and B is the identity matrix. The 
matrix B has been included in this investigation for completeness. 

The approximation to the quadratic cost function in Eq. (3.2) becomes 

(3.11) J N (0$,& N ) = J o + <g W (')^ ,V g ,V (0)]-* 

where 0^ and R N are now matrices which weight the state vector oo^ and the control input g w respectively. 

The finite dimensional control problem is stated as 



subject to the governing equations in Eq. (3.10). The optimal control input which minimizes Eq. (3.12) is given by 


(3.13) Jw/wM = -(R N )~ l B N *n N u> N (t) 

where fl^ is the solution to the finite dimensional ARE 


(3.14) 


a n * n* + n N A N - n N B N (R N y 1 B N *n N + q n = o 


As in the infinite dimensional case, a solution F/ exists for Eq. (3.14) it the pair (A a ,B n ) is stabilizable and the pair 
{A N ,C N ) is detectable [4]. This is easily satisfied since actuation and sensing are performed everywhere in the flow 
field. 

Equation Eq. (3. 1 3) shows that the finite dimensional approximation to the optimal feedback gain K opt is 
(3.15) Kop = -{R n )~ { B n *U n 

Note that the finite dimensional approximation K N to the infinite dimensional feedback control K is dependent on the 
approximation to the ARE solution F/ and the adjoint B N \ In this study, the B w matrix is /jvxN which means that the 
approximation K N is primarily dependent on the approximation P1 A . 
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4. Reduced Order Model. The computational cost of solving Eq. (3.14) is on the order of N 3 floating-point 
operations (flops), where N is the number of states. For two dimensional flows the number of flow states, for a 
computational simulation, can range anywhere from 10 4 to 1(>\ For a simple two dimensional flow problem, say 
N = I0 4 , the cost would be on the order of I0 12 flops. For a larger two dimensional flow problem, say N = l(>\ 
the computational cost would increase by a factor of 1000, to 10 15 flops. For three dimensional flows the number of 
states can range from 10 5 to 10 6 . Thus, the cost of solving the ARE can grow prohibitively expensive, even for two 
dimensional flow problems. 

The computational cost can be reduced significantly by using a Chandrasekhar system approach |3]. The only 
draw back to the Chandrasekhar system is that the number of inputs must be much smaller than the number of states 
in order to significantly reduce the cost of solving the Riccati equation. Since the number of inputs for this problem 
are equal to the number of states, a different approach needs to be taken. 

Another way to reduce the cost, and the approach taken here, is to develop a reduced order model of the system 
given in Eq. (3.10). This reduced order model has a smaller number of states which reduces the cost of solving the 
Riccati equation in the LQR control design. Since the goal is to approximate the infinite dimensional LQR feedback 
gain, K opt , a reduced order model which gives a good approximation should be chosen. It would also be desirable for 
the approximation to the infinite dimensional gain to converge in some reasonable way as the number of states in the 
reduced model increase. 


4.1. Krylov Space Method. Most model reduction methods for feedback control systems use a Hankel-norm 
approach. This approach calculates Hankel singular values based on the controllability and observability of a given 
system. The advantage of this approach is that the norm of the modeling error is bounded by the sum of the Hankel 
singular values not retained in the reduced order model. The draw back it that the cost of this approach is of order N*. 
In this study a Krylov based method is used to project Eq. (3.10) onto a reduced Krylov subspace which includes the 
leading modes of the larger system. This approach is significantly cheaper than the Hankel-norm approach but does not 
take into account the controllability and observability of the closed-loop system. It can be argued, that since actuation 
and sensing is performed everywhere in the flow field, that the leading modes of the homogeneous system are the 
dominant modes for the closed-loop system. If the matrix B N had some structure other than identity and sensing was 
not performed everywhere, then the leading modes of the homogeneous system are not necessarily the most dominant 
modes (modes with the largest real part) for the closed- loop system. 

Consider the homogeneous solution to Eq. (3. 10), 


(4.1) vf i {t) = T{t)<s% 

. l\ f 

where T(t) = e* ' is a Co-semigroup. Equation (4. 1 ) can also be written as. 


(4.2) 


co" (kAt + At) = T(kAt + A/)o# 
= T(At)T(kAt)v$ 
= T(At)o/ i (kAt) 
o/ +l = T(At)(rf 


where (O* = CO' v (kAt) and the superscripts N have been dropped for convenience. Using a Krylov method, the action 
of the semigroup operator T(At) can be approximated by a reduced Krylov subspace. A /.'''-dimensional Krylov 
subspace, given an operator T(At) and a vector v, is defined as 


(4.3) 


3u(7\v) = s/?a/i{v,7'v,7' 2 v,...,7'* 1 v} 
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This Krylov subspace is generated using the ARPACK software package which uses an Arnold i/Lanczos scheme with 
implicit restarts [19]. The orthonormal basis for the Krylov subspace is generated by following the Arnold! Process [ 1 ]. 

V; = W/Vllw/H 

(4.4) k 

Wy+l = Twj - £ \ m (v m ,T\j) 

m= 1 

The vectors vj for j = 1, 2,..., A are computed using a given starting vector w\. The resulting NxK vector V = 

( v i , vt is known as a Ritz Vector where V € 3Q(7>, ). Note that the generation of the Krylov subspace only 
requires the action of T on a given vector and not the explicit matrix T. As shown in Eq. (4.2) the action of T on a 
given vector is just the advancement of the vector, or in this case the vorticity field, by one time step A/. 

Equation (4.4) can be restated in the standard matrix form of the Arnoldi decomposition as 

(4.5) TV = VH + w* + ie[ 

where w*+ie[ is a remainder term. The matrix H is a Kx K upper Hessenberg matrix. The operator VV T is a Ax A 
projection operator onto the Krylov subspace and V 1 V is the K x K identity matrix. 

The action of the semigroup operator T can be approximated by projecting it onto a Krylov subspace 

(4.6) T ss VV t TVV t « VHV 7 
where H « V t TV for small w^+ie^ . Diagonalizing the matrix H yields 

(4.7) H = EAE~' 

where A is a diagonal matrix of eigenvalues and £ is a matrix whose columns are the corresponding eigenvectors. 
The Krylov subspace is generated by using implicit restarts which keep the leading eigenvalues of the T semigroup 
operator. Using the definition of 7" (A/) = ^ and Eq. (4.6) an approximation to the operator/) results in 

TiAt)-^ « VHV 1 

(4.8) 

= VEAE-'V 1 

Solving for A produces 

(4.9) 

where H r — E log(\)E~* / At. 

4.2. Control Problem. Projection of the A"' dimensional system in Eq. (3.10) onto a k"' dimensional Krylov 
subspace results in the reduced system 

z = H r z + B r g N 

(4.10) 

z(0) = z 0 

where B r = V T and the aggregated state vector z is defined as 

(4.11) z = V T (a N 
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Thus, the state vector 0)^ has been approximated by another state vector z, constrained to stay in the Krylov subspace 
given by the basis V 7 . The vorticity state vector is then approximated as 0) % Vz. 

The system in Eq. (4. 10) is now used to compute an approximation to the the desired optimal feedback gain K N 
given in Eq. (3.13). Using the Krylov subspace V the cost function in Eq. (3.1 1) results in 

(4.12) J(z 0'tP) = f Q [{z(t),Q r z(t)) + {g N (t),Rrg N (t))]dl 

where R r = and Q r — V r Q N V. The control problem defined in Eq. (3.12) for the reduced system becomes 

(4.13) min7(zo,g jV ) 

g v 

subject to the governing equations in Eq. (4. 10). The optimal control input for this problem is 

(4.14) g ? ll , l (t) = -R; l B;n r z(t) 
where the (k x k) matrix, n r , is the solution to the finite dimensional ARE 

(4. is) h; n r + n r H r - n r B r R; 1 + q, = o 

Using Eq. (4. 1 1 ) and Eq. (4. 1 3), the approximation to the desired feedback gain K N becomes 
(4.16) K n « ~R~ ] vn,v r 


where the adjoint B' r = V. 

5. Numerical Method. The calculation of T(At) on a given vector is achieved by computing a time accurate 
solution to the linearized Navier-Stokes equations described in Eq. (2.8). A solution to Eq. (2.8) is computed using a 
semi-implicit finite difference scheme. This scheme approximates the spatial derivatives using a second-order central 
difference method for the viscous term and a third-order upwind scheme for the convection terms. These equations 
are then solved using a multigrid acceleration method with Gauss-Seidel relaxation. 


5.1. Time Discretization. The governing equations in Eq. (2.8) are discretized in time using an explicit up- 
winding scheme on the convection terms and a implicit Crank-Nicolson scheme on the diffusive term. This scheme is 
locally second-order accurate in space and first-order in time. The discretization of Eq. (2.8) has the form 


(0 


,/j+i 


Ki 


-V" - (Qu")„ - V" • (co"U),, + ~v 2h (<t' + (o;: y ) 


(5.1) 


At 

v 2/ V;:, 


= — (D 




where to)' • = (ti(ih,jh,nAt). The operators V ; '. and V 2 ^ are discrete approximations to the operators V and V 2 , respec- 
tively. The discretized vorticity transport equations in Eq. (5. 1 ) can now be expressed as 


(5.2) 




oo 




= (0 lj - v h ■ (Qu n )/J - V" ■ (co"U )ij + 


2Re 


V 2, 'to'’ 


The nonlinear Navier-Stokes equations in Eq. (2.1) are similarly discretized by replacing D. with of' and U with u n 
in the equation above. The steady state solution to the nonlinear equations in Eq. (2.1 ) are used in the linearized 
Navier-Stokes equation for the base flow field. 

The discretized equation in Eq. (5.2) can be expressed in the form 


(5.3) 


Ww" +I =b 
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where b is the right hand side of Eq. (5.2) and M the discrete operator (/ — At/2Re V 2/, )> Since M is a large and 
sparse matrix, the solution for co" +1 favors an iterative method. To accelerate the iterative method, a multigrid routine 
which uses Gauss-Seidel relaxation is used. Likewise the solution to the Poisson equation in Eq. (5.1 ) is solved using 
an iterative method with multigrid acceleration. 

5.2. Spatial Discretization. The Laplacian operator V 2 is discretized using a standard second-order central dif- 
ference scheme and has the form 


(5.4) 


V" 


->h 


“5 = 


co;y IJ -2(o;; / + < 


1,7 


CO" - 2(0" • + CO" 


Av J 


+ 


'■/ + 1 




u >,i - 1 


Ay 2 


The first order derivatives in the linear convection terms are evaluated using a four-point upwind scheme. The x 
derivative component for the term V h • (ftu n );j, for u > 0, is approximated by the upwind scheme 


(5.5) 


d(Qu)ij _ -(fl«) ;+ | J + 3(fl//) lJ - 3(fl«),_| / + 2(Qu) i+iJ - (flu),-. 


d.x 


3Av 


1,7 


2Av 


+ 0( Ar 


where the parameter q controls degree of modification to the central difference term [6]. If q — 0.5 then the scheme 
in Eq. (5.5) becomes an upwind scheme of 0(Ar 3 ). If q = 0 then the scheme is reduced to a second-order central 
difference scheme. The central difference approximation for the linear convection term has good accuracy but will 
produce oscillations when the mesh Peclet number (p = u dx Re) becomes greater than 2. The upwind scheme will 
reduce these oscillations for p > 2 but at a cost of reduced accuracy. This reduced accuracy is a result of artificial 
diffusivity added by the upwind scheme. Therefore, a switch is used which sets q — 0 when p, . j < 2 and q — 0.5 when 
p, j > 2. This switch results in a scheme which has better accuracy at low mesh Peclet numbers than a straight upwind 
scheme, yet retains the advantages of an upwind scheme for large Mesh Peclet numbers. 

5.3. Boundary Conditions. The boundary condition for the no-slip walls are computed using Jensen's formula- 
tion [9] attributed to Jensen by Roache [ 17]. This formulation, also known as Briley’s formulation and was used by 
Pearson [14] and Ghia et al. [7]. 

Jensen’s formulation computes a boundary value for the vorticity by taking a Taylor series expansion of the stream 
function normal to the wall. Jensen’s formula is given by 


(5.6) 


oooj = 


7\|/o,;-8y]j + i|/2j 

2 h 2 


3ho,j 

h 


+ o(C~) 


This boundary condition is claimed by [9, 17, 14, 7] to be 0(/? 2 ) for the vorticity at the boundary. It has been shown 
in [20] that this formulation can be thought of as an 0(/r 3 ) approximation to the Neuman boundary condition on the 
stream function in Eq. (2.9b) rather than a 0(h 2 ) boundary condition for vorticity. If Eq. (5.6) is rewritten as 


(5.7) 


h 

3«o ,/ 


7 — &Vij + \p2j 

6 h 


— wo + O (/r 3 ) 


then the lim/,^o recovers the equation (3^/5«)|o,/ = —mo,/. This is consistent with Eq. (2.9b), which is the equation 
being modeled. Therefore the condition on the vorticity at the boundary can be thought of as an 0{hr) approximation 
to Eq. (2.9b) instead of an 0(h 2 ) boundary condition for vorticity. Spotz used Jensen’s formula using a compact 4 th 
order method and showed that this formula resulted in an 0(/r 3 ) approximation. 


5.4. Driven Cavity Problem. A numerical simulation of the driven cavity problem, using the full nonlinear 
Navier-Stokes equations in Eq. (2.1 ), is used to give a measure of validation for the proposed numerical method. The 
driven cavity problem is a typical two-dimensional model problem that is used to evaluate and compare numerical 
methods for incompressible viscous flows. Most notable are the steady state results published by Ghiaet al. [7]. 
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Ik i 5.2. Stream lines of the driven cavity problem at steady 
state for Re= 1 000, 


FlG. 5.3. A contour plot of the steady state vorticity- field for 
the driven cavity problem at Re -J 000. 


The driven cavity problem, shown in Fig. 5.1, has a top wall which moves at a nondimensional velocity of 
Ujop = 1 • The moving wall induces the flow as a result of the viscous forces generated at the moving wall. The steady 
state solution is then computed on a uniform grid which has 129 grid points in both the x and v directions. A steady 
state flow field is found by simply marching the impulsively started cavity in time, until a satisfactory steady state 
solution is achieved. Its interesting to note that this problem has two singularities at the top two corners where the top 
moving wall meets the two stationary side walls. The impact from these singularities are considered to be small since 
the grid is relatively fine. 

Figure 5.2 and 5.3 show the stream lines and vorticity contours for the steady state solution at Re = 1000. The 
stream lines in Fig. 5.2 show two recirculation zones at the bottom comers of the cavity. The size and location of these 
recirculation zones compare very well the numerical results given by Ghia et al. [7]. The minimum stream function 
was computed to be ^ min — 0.1 17985 and compared very well to the minimum stream function, ^ min = 0.1 17929, 
computed by Ghia. The location of the center of the main vortex is computed to be at x = (0.5312,0.5625) as 
compared to x = (0.5313.0.5625) reported by Ghia. The contour plot of the vorticity field in Fig. 5.3 also compares 
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Fig. 5.4. A comparison of the horizontal velocity profile, v, 
at the centerline of the carin', \ = 0,5. with the profiles reported by 
Ghia et al. \7j at Re= 100. 400. I (XX), 32(H). 


Fig. 5.5. A comparison of the vertical velocity profile, u. at the 
centerline of the cavity. \ — 0.5 with the profiles reported by Ghia et 
al. 17J at Re = 100,400, 1000,3200. 


very well with the results from Ghia. 

Figure 5.4 shows the horizontal velocity profiles, v\ at v — 0.5 for Re — 100,400, 1000, and 3200. These velocity 
profiles are compared to the velocity profiles reported by Ghia et al. [ 7 J . Likewise the vertical velocity profiles for 
u at A' — 0.5 are shown in Fig. 5.5 and then compared to the numerical data given by Ghia. These two figures show 
good agreement between the current results and the data from the Ghia paper. Note that a 129x 129 grid was used for 
all the Reynolds number cases except for the Re — 3200 case where a 257 x 257 grid was used as was done by Ghia. 
Based on these results a measure of validation has been shown for the proposed numerical method for the calculation 
of steady state flows. 

6. Results. Application of the distributed feedback control design is now applied to an unstable shear layer 
problem and a stable driven cavity problem. The computed eigenvalues and eigenvectors for the bounded shear layer 
flow will be compared to eigenvalues and eigenvectors from an Orr-Sommerfeld analysis. This comparison will 
provide a measure of validation for the reduced order model. 

6.1. Shear Flow Problem. The method described above is now applied to a two-dimensional bounded shear 
layer problem as shown in Fig. 6.1. In this problem there are two layers of parallel fluid traveling in opposite directions. 
At the intersection of these two flows is a shear layer which has a hyperbolic tangent velocity profile. The large velocity 
gradient in the shear layer results in a large concentration of vorticity. This type of flow pattern is inviscidly unstable 
to small disturbances. The base flow field for this problem is described by the equations 
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Fl(i. 6. 1 . The velocity' profile in the v direction for the shear layer problem. 


( 6 . 1 ) 


f/(.v,v) = [l/(y),0] = [f/otanh(y//>),0] 


Q(jc,v) 


sech 2 (v/fc) 

b 


(-1 <*< 1,-1 < y < 1 ) 


where /? = 1 / 30, Re = 100 and the flow is periodic in the x direction. The Reynolds number for the shear layer 
problem is typically defined using the characteristic length /?, which determines the width of the shear layer. The 
Reynolds number for the shear layer is defined as 


( 6 . 2 ) 


Re = 


Unb 


v 


where Uo is the nondimensional characteristic velocity and v the kinematic viscosity of the fluid. The U(y) profile is 
shown in Fig. 6. 1 where the boundaries moving walls with a no-slip boundary condition. 

The velocity profile for this problem represents an exact solution to the steady state inviscid flow equations. 
However it does not exactly satisfy the steady state equations for the incompressible viscous flow. The unsteadiness, 
in this base flow, is due to the viscous diffusion term which wants to diffuse the high concentration of vorticity in 
the shear layer. This causes the vorticity layer to spread out and diffuse away from the center region. In practice this 
viscous diffusion effect is overlooked in the formulation of the linear stability problems. [15] 


6.1.1. Reduced Order Model. Using the shear layer base flow field given in Eq. (6.1), a Krylov basis vector 
V of dimension K = 201 was constructed to form the reduced order model in Eq. (4.10). These vectors were found 
using the ARPACK software package where the top 201 eigenvectors with the largest real part were computed. This 
program starts by generating a large Krylov subspace of 600 Ritz vectors. It then performs an implicit restart which 
removes unwanted Ritz vectors and then generates new Ritz vectors, replacing the discarded vectors. This process 
is repeated until the desired number of leading eigenvectors has converged. For this problem, 2735 T 0 ) operations 
(time steps) were computed by the flow solver. The finite difference approximation was made using 64 grid points 
in the streamwise direction and 128 grid points in the cross stream direction. Figure 6.2 shows the convergence of 
the ARPACK routine where an implicit restart was performed at the start of a new iteration. The convergence rate is 
dependent on the number of desired eigenvectors, the size of the Krylov subspace chosen, and the size of the time step. 

The top 201 eigenvalues, with the largest real parts, are shown in Fig. 6.3. These eigenvalues are compared to the 
eigenvalues computed from an Orr-Sommerfeld stability analysis using the spectral method presented by Orszag[12]. 
The numerical method used by Orszag was derived for the stability of plane Poiseuille flow which was then modified 
for the shear layer problem. The comparison between the eigenvalues shows how the Krylov method was able to corn- 
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Fi(i, 6.2. The convergence history of the Krylov vectors usin g the AR PACK software package. The Krylov subspace had 600 vectors and was 
looking for the top 201 eigenvectors with the largest real parts . 

pute the eigenvalues of the shear layer problem reasonably well. From this figure it can be seen, that the eigenvalues 
start to degrade as the imaginary parts become larger and as the real parts of the eigenvalues become increasingly 
negative. The reduced accuracy of the eigenvalues in these areas show the limitations of the finite difference approxi- 
mation and the computation of the eigenvalues using the Krylov time stepping method. As the imaginary and negative 
real parts of the eigenvalues increase, so do the spatial oscillations of the eigenvectors. Therefore the spatial resolution 
of the finite difference method governs how many of the eigenvectors can be resolved accurately. For the shear layer 
problem, the number of modes that can be resolve are sufficient our investigation. 

A sample of the eigenvectors computed using the Krylov method and the Orr-Sommerfeld analysis are shown in 
Fig. 6.4. A Fast Fourier Transform (FFT) was taken of the eigenvectors, generated by the Krylov method. This FFT 
showed that the eigenvectors have one dominant wave number in the x direction. Therefore the eigenvectors can be 
expressed in the form 

(6.3) V(x ,y) = v(y)e iox 

where Fig. 6.4 is showing the function v(y) for four different eigenvectors. Figure 6.4 a shows the most unstable 
eigenvector for the shear flow problem. This mode has a wave number of (5 — In and shows the form of the instability 
of the shear layer. This eigenvector compares very well to the eigenvector computed using the Orr-Sommerfeld spectral 
analysis. This unstable mode is the easiest mode to capture using this Krylov method since it grows very fast. The 
other eigenvectors in Fig. 6.4/? through d are stable modes. The eigenvectors shown in Fig. 6.4/? and d both have 
wave numbers o = 7t and the eigenvector in Fig. 6 Ad a wave number of o = 2ji. These modes compare very well with 
the corresponding eigenvector from the Orr-Sommerfeld analysis. This comparison shows a measure of validation 
in computing the leading eigenvectors and eigenvalues using the Krylov method. This comparison also shows that 
the eigenvalues start to degrade in accuracy for eigenvalues with increasing wave numbers and negative real parts. 
The accuracy of the eigenvalues and eigenvectors can be improved by increasing the spatial resolution and by using 
a numerical scheme of higher order. It was also seen that decreasing the time step could improve the accuracy of the 
eigenvalues up to some limit. Further reduction of the time step beyond this limit did not increase the accuracy of the 
eigenvalues. 

6.1.2. LQR Control Design. The optimal feedback kernel for the reduced system was computed by minimizing 
the quadratic cost function in Eq. (4.12). The weights for the cost function in Eq. (3.1 1 ) were set to R = / and Q = /. 
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This resulted in the weights for Eq. (4. 1 2) 

(6.4) R r = /, Q r = V t QV = V T V - 1 

The ARE for the reduced system, given by Eq. (4.15), is then solved for n r using the MATLAB software package. 
This reduced system results in an ARE solution with 201 states, which is large for typical ARE problems, but is 
significantly smaller than the 8192 states of the full system. Figure 6.5 shows the diagonal elements of the ARE 
solution U r as a function of the mode number. The modes are the eigenvectors of the linearized flow equations and the 
eigenvalue Xj corresponds to the i th mode where Re(Xj) > Re(Xj+ 1 ). This figure shows how the first six modes, which 
are unstable modes, result in the largest contribution to the feedback kernel K N and how the contribution decreases for 
the higher mode numbers. This figure basically shows the convergence of the feedback kernel with respect to the size 
of the reduced order system. 

Substituting n r into Eq. (4.16) results in the approximation to the finite dimensional feedback kernel K N . Using 
the approximation to K N the optimal feedback control can be computed as 

(6.5) g»,„(0 = K N aP{t) 

where K N is a (8256x8256) matrix, g„ pt a (8256 x 1) column vector, and a (1 x8256) row vector. It can be 
seen from Eq. (6.5) that the n lh column of K N corresponds to the distributed control input g for a disturbance to = e„. 
Similarly, the m th row of K N correspond to the distributed feedback gain for the control input at g m . To illustrate 
this idea Fig. 6.6 shows the distributed control ; given a unit disturbance at four different locations in the flow. 
Figure 6.6 a shows the point disturbance on the bottom wall and shows a local positive control near the point distur- 
bance. Fig. 6.6 b shows the disturbance just off the wall at y — —0.75 and shows the localized nature of this feedback 
control. In Fig. 6.6c the disturbance is move closer to the shear layer. This figure shows the same type of localized 
feedback control around the disturbance as before but with a large control force. The disturbance near the shear layer 
also shows an interesting feedback control force which appears in the shear layer. The point disturbance is then moved 
to the center of the shear layer as is shown in Fig. 6.6 d. This distributed feedback control shows a pattern which is 
similar to the unstable modes in the flow. This should be expected since the unstable modes were associated with the 
largest feedback gains computed in ll r . 

In an effort to quantify the spatial structure of the control effort a norm is defined in the following way 

X (O 2 

m~ 1 

where c„ is the n th element of a vector defining the control effort and c — (ci ...,cw). The idea is for c n to be 
a measure of the control effort in an Li norm sense for a point disturbance located at co = e„. Likewise the spatial 
structure of the feedback gain can be evaluated by defining the norm 

1 1/2 

(O 2 

where f — (/i ,/is ...,/v) is a vector describing the measure of the feedback effort. The value of /,„ is an Lj norm of 
the distributed feedback gain for actuation at g, n . Note that the feedback kernel is symmetric since the weight R and 
the matrix B are symmetric. Therefore the the vectors c and f, which describe the control effort and feedback effort, 
are equal. 

A plot of the control effort at a constant x value, for the shear flow problem, is shown in Fig. 6.7. This figure 
shows how the control effort, and consequently the feedback gain, is large in the shear region. This figure also shows 
that the shear region would be most favorable for actuation and sensing as might be expected. 
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6.2. Cavity Flow Problem. The second problem involved the design of a distributed feedback controller for a 
driven cavity problem. Unlike the shear layer, the base flow field for the driven cavity was stable and satisfies the 
finite difference approximation to the steady state flow equations. Since the flow is stable, the model will not contain 
any unstable modes as was seen in the shear layer problem. This means that the cavity problem will not have a small 
number of dominant modes in the reduced order model and feedback kernel as was seen in the shear layer problem. 

6.2.1. Reduced Order Model. The base flow field for the cavity problem was the steady state flow field com- 
puted in section 5.4 for Re = 1000. The vorticity for the base flow is shown in Fig. 5.3 and the stream function in 
Fig. 5.2. Using this base flow field, a Krylov basis vector V of dimension K = 400 was constructed to form the reduced 
order model of the linearized flow equations. Figure 6.8 show the convergence of the top eigenvalues with the largest 
real parts. It took 50 implicit restarts and 15370 T 0) operations (time steps) by the flow solver in order for the top 400 
modes to converge. The distribution of the top eigenvalues is shown in Fig. 6.9. A sample of the eigenvectors with the 
largest real parts are shown in Fig. 6.10. This figure shows the real and imaginary parts of the eigenvectors. Notice 
that some of the eigenvectors are real and do not have any imaginary parts. 

6.2.2. LQR Control Design. Using the reduced order system for the linearized cavity flow, the optimal feedback 
kernel is then computed using an LQR control design approach. As in the shear layer problem, the Rieeati solution 
U r to the ARE in Eq. (4.15) was computed which minimizes the cost function in Eq. (4. 12). The weights for the cost 
function in Eq. (3.11) are uniform and set to R = / and Q = / which results in the weights R, = / and (?, =/. Figure 
6.1 1 shows the diagonal elements lor n r as a function of the mode number. This figure shows a similar decay rate for 
the diagonal elements of the Riccati solution H, as compared to the shear flow problem. 

The feedback kernel K N can now be approximated by substituting the Riccati solution to the reduced system 
n r into Eq. (4.16). This results in a ( 1664 1 x 16641 ) feedback kernel K N where the approximation to the optimal 
distributed control g^ p , is given by 

(6.8) g = K N wf*(t) 

Using this optimal feedback kernel, the distributed control for a point disturbance can be computed. Figure 6. 1 2 shows 
the approximation to the distributed control given a point disturbance at four different locations in the cavity. A point 
disturbance on the bottom wall is shown in Fig 6. 1 2 a with the resulting distributed control force. Figures 6. 1 2b through 
d show the point disturbance in the interior of the cavity away from the effects of the wall. The distributed feedback 
for these interior points show a smooth circular pattern around the disturbance. This figure also shows that the applied 
control force is a maximum at the point disturbance and then decays as it moves away from the disturbance location. 
It is interesting to see how the control force decays as it moves away from the location of the point disturbance. 
Figure 6.13a shows a mesh plot of the distributed control shown in Fig 6.12r and Fig 6.13/? shows a plot of the 
distributed control for y = 0.5. These figures illustrate how the control force decays from the point disturbance located 
at (.v.v) = (0.5, 0.5). 

The distributed control effort and feedback gain for the cavity problem is shown as a mesh plot in Fig. 6.14 
and as a contour plot in Fig. 6.15. These two figures show how the distributed control effort and feedback gain are 
concentrated in the center of the cavity. There are also some peaks near the center of the walls and one large peak 
where the flow induced from the top moving lid impinges on the right stationary wall. The figure also shows that the 
control effort is smaller at the corner of the cavity. Therefore sensing and actuation would be most favorable near the 
center of the cavity, near the center of the wall, and on the top part of the right wall. 

7. Conclusion. This study has demonstrated how a Krylov subspace method can be used to derive a reduced 
order model of the linearized incompressible Navier-Stokes equations and applied to a two-dimensional shear flow 
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and driven cavity problem. By assuming sensing and actuation everywhere in the flow field an optimal feedback 
kernel can be found. This feedback kernel provides insight into the regions of the flow where the control effort is 
Marge’ and where the feedback gain is Marge’. This information can be used to isolate the regions of the flow field 
where sensing and actuation are most favorable. This information also shows where one should be searching for the 
best place to put actuators and sensor. This knowledge has the potential to reduce the search areas of the flow field for 
the placement of actuators and sensors. In the shear layer case, the search area could be significantly reduced, where 
as the driven cavity problem only showed a slight reduction in the search area was achieved. 

The two-dimensional shear flow problem showed that sensing and actuation was most favorable in the shear layer 
as might be expected. The driven cavity problem showed that a region in the center and parts of the walls were the 
most favorable tor the placement of sensors and actuators. It also showed that there was not a dominant region, as 
in the shear layer problem, and that the corners of the cavity and a region near the walls were unfavorable for the 
placement of sensors and actuators. 

The advantage of this approach is that a simple time stepping vorticity stream function code could be used to 
derive the linearized model of the incompressible Navier-Stokes equations. To improve on this approach, higher-order 
spatial discretizations should be used to improve the spatial accuracy of the higher frequency modes. The drawback 
to the time stepping approach used here is the convergence of the eigenvalues and eigenvectors. To solve this problem 
a shift invert approach described by Sorensen [19] can be used but this requires an explicit representation of the flow 
equations as oppose to the coupled vorticity stream function equations used here. 
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Mg . 6. V A comparison of the 201 eigenvalues computed using the Krylov method using the two-dimensional finite difference numerical 
approximation to the one-dimensional Orr-Sommerfeld analysis using a spectral method at Re = 100. The eigenvalues for the Orr-Sommerfeld 
equations were computed using wave numbers . in the x direction, equal to 0, k,2kJtz, and 4k. 



Pig. 6.4. A comparison of eigenvectors computed using the Krylov method to a one-dimensional Orr-Sommerfeld analysis which is solved 
using a spectral method. 
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Fig. 6.5. The value of the diagonal elements for U r . the solution of the ARE for the reduced system. The first six modes are unstable and 
therefore produce large values in n, solution. The other modes are all stable and the contribution of these modes to the feedback kernel decreases 
for the modes with increasing negative real parts. 



Fig. 6.6. This figure shows the spatial distribution of the applied control ^ v (x) for a unit disturbance which is shown by the black dot. In 
these four figures the unit disturbance is moved along the v axis from the w all to the center of the shear layer. 
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Pic*. 6.7. A measure of the control effort and feedback gain. 



Fici. 6.8. The convergence history' of the Krylov vectors using the ARPACK software package for the cavity ’ problem. The Krylov subspace 
had 8(X) vectors and was looking for the top 400 eigenvectors with the largest real parts. 
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Pic. 6.9. Eigenvalues of the driven cavin' problem for Re — 1000. 
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Fig. 6 . 1 5. A contour plot showing the control and feedback effort for the forced cavity problem. 
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